R Notebook

library(stringr)
library(tidyverse)
library(bigrquery)

response <- try(system('~/google-cloud-sdk/bin/gcloud projects list --quiet', intern = T))
projectid <- strsplit(response[2], " ")[[1]][1]

options(na.action = "na.fail") 

source("./dredge_functions.R")

dredge_and_subset <- function(data) {
  model <- lm(
    presence_ratio ~ foraging_niche + trophic_niche + is_noctural + pc1 + pc2 + pc3 + pc4 + iucn_red_list_2020 + realm,
    data=data
  )
  dredge_result <- dredge(model)
  summary(model.avg(dredge_result, subset = delta < 10))
}

load_species_data <- function(poolname) {
  filename <- str_replace('species_metrics_##POOL_NAME##.csv', '##POOL_NAME##', poolname)
  
  if (!file.exists(filename)) {
    column_name = str_replace_all('present_in_##POOL_NAME##_pool', '##POOL_NAME##', poolname)
    if (poolname == 'both') {
      column_name = 'present_in_both_pools'
    }
    
    sql <- str_replace_all("
         WITH species AS (
              SELECT 
                  scientific_name,
                  count(*) AS regional_pool_count,
                  countif(present_in_city = true) AS city_count
              FROM model.regional_species_filtered
              JOIN model.taxonomy USING (scientific_name)
              WHERE ##COL_SELECT## = true
              GROUP BY scientific_name
          )
          SELECT 
              species.*,
              body_morphspace.pc1,
              body_morphspace.pc2,
              body_morphspace.pc3,
              body_morphspace.pc4,
              trophic_niche,
              foraging_niche,
              is_noctural,
              iucn_red_list_2020,
              realm,
              ttc.cluster,
              ttc.cluster_half,
              ttc.cluster_double
          FROM species
          JOIN model.taxonomy USING (scientific_name)
          JOIN datasci.taxonomic_trait_clusters ttc USING (scientific_name)", '##COL_SELECT##', column_name)
  
    tb <- bq_project_query(projectid, sql)

    data <- bq_table_download(tb)
    write_csv(data, filename)
  }
  
  data <- read_csv(filename)
  
  data[is.na(data$foraging_niche),]$foraging_niche <- 'Omnivore'
  
  data$foraging_niche = as.factor(data$foraging_niche)
  data$trophic_niche = as.factor(data$trophic_niche)
  data$is_noctural = as.factor(data$is_noctural)
  data$presence_ratio = data$city_count / data$regional_pool_count
  data$iucn_red_list_2020 = as.factor(data$iucn_red_list_2020)
  data$realm = as.factor(data$realm)
  
  data
}

merlin_data <- load_species_data('merlin')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  iucn_red_list_2020 = col_character(),
  realm = col_character(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

merlin_data

merlin_result <- dredge_and_subset(merlin_data)

Fixed term is "(Intercept)"

merlin_summ <- model_summary('merlin', 'species', merlin_result)
merlin_summ

birdlife_data <- load_species_data('birdlife')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  iucn_red_list_2020 = col_character(),
  realm = col_character(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

birdlife_result <- dredge_and_subset(birdlife_data)

Fixed term is "(Intercept)"

birdlife_summ <- model_summary('birdlife', 'species', birdlife_result)
birdlife_summ

both_data <- load_species_data('both')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  iucn_red_list_2020 = col_character(),
  realm = col_character(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

both_result <- dredge_and_subset(both_data)

Fixed term is "(Intercept)"

both_summ <- model_summary('both', 'species', both_result)
both_summ

either_data <- load_species_data('either')

── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  scientific_name = col_character(),
  regional_pool_count = col_double(),
  city_count = col_double(),
  pc1 = col_double(),
  pc2 = col_double(),
  pc3 = col_double(),
  pc4 = col_double(),
  trophic_niche = col_character(),
  foraging_niche = col_character(),
  is_noctural = col_logical(),
  iucn_red_list_2020 = col_character(),
  realm = col_character(),
  cluster = col_double(),
  cluster_half = col_double(),
  cluster_double = col_double()
)

either_result <- dredge_and_subset(either_data)

Fixed term is "(Intercept)"

either_summ <- model_summary('either', 'species', either_result)
either_summ

Full result

all_species_results <- full_join(full_join(merlin_summ, birdlife_summ), full_join(both_summ, either_summ))

Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")

write_csv(all_species_results, "species_result.csv")

merlin_data$source <- 'Merlin'
birdlife_data$source <- 'Birdlife'
both_data$source <- 'Both'
either_data$source <- 'Either'

plot_data <- rbind(merlin_data, birdlife_data, both_data, either_data)

ggplot(plot_data, aes(x = pc1, y = pc2, alpha = presence_ratio, color = as.factor(foraging_niche))) + geom_point() + theme(legend.position = "none") + facet_wrap(~ source)

ggplot(plot_data, aes(x = pc3, y = pc4, alpha = presence_ratio, color = as.factor(foraging_niche))) + geom_point() + theme(legend.position = "none") + facet_wrap(~ source)

library(ggpubr)

merlin_data$niche <- paste(merlin_data$trophic_niche, ' ', merlin_data$foraging_niche, ' ', merlin_data$is_noctural, ' ', merlin_data$cluster)

gg <- function(gg) {
  gg + geom_point() + theme_bw() + labs(alpha="Presence Ratio", colour="Niche")
}

myplot <- ggarrange(
  ncol = 2,
  nrow = 2,
  gg(ggplot(merlin_data, aes(x = pc1, y = pc3, alpha = presence_ratio, colour = niche)) + ylab("Short Tail + Pointy Beak -> Long Tail + Stubby Beak")) + theme(legend.position="none", axis.title.x=element_blank()) ,
  gg(ggplot(merlin_data, aes(x = pc2, y = pc3, alpha = presence_ratio, colour = niche))) + theme(legend.position="none", axis.title.x=element_blank(), axis.title.y=element_blank()),
  gg(ggplot(merlin_data, aes(x = pc1, y = pc4, alpha = presence_ratio, colour = niche)) + ylab("Long Tail + Pointy Beak -> Short Tail + Stubby Beak") + xlab("Body Size")) + theme(legend.position="none"),
  gg(ggplot(merlin_data, aes(x = pc2, y = pc4, alpha = presence_ratio, colour = niche)) + xlab("Beak Size"))  + theme(legend.position="none", axis.title.y=element_blank())
)
myplot <- annotate_figure(myplot, top = text_grob("Species Presence based on Morphology", face = "bold", size = 14),
                          bottom = text_grob("Colour: Niche, Alpha: Presence Ratio", hjust = 1, size = 10))
myplot

jpeg("species.jpg", width = 1024, height = 1024)
myplot
dev.off()

null device 
          1